GS <- mutate(GS,Date = as.Date(Date, format = "%d-%b-%y"))
GSxts <- tk_xts(GS)
## Warning in tk_xts_.data.frame(data = data, select = select, date_var =
## date_var, : Non-numeric columns being dropped: Date
## Using column `Date` for date_var.
#pick the earliest record date as the start date of the single trade
startD <- min(index(GSxts))
#end date of the single trade is after 30 calendar days
endD <- startD+30
#adjust the end date backwards if end date (a calendar day) is not in the xts
while(!endD %in% index(GSxts))
  endD <- endD-1

xts_obj <- GSxts[paste(c(startD,endD),collapse = "/")]
quantity = 100
dates <- index(xts_obj)
start_date <- min(dates)
end_date <- max(dates)
start_price <- as.numeric(xts_obj[start_date, "Close"])
start_volatility <- as.numeric(xts_obj[start_date, "IV30"])

df <- tibble(Date = dates)
df$Close <- coredata(xts_obj[, "Close"])
df$IV30 <- coredata(xts_obj[, "IV30"])
r <- 0.8 / 100
X <- start_price/(exp(qnorm(0.25)*start_volatility/100*sqrt(30/365) - (r+0.5*(start_volatility/100)^2)*30/365))
sigma = start_volatility

# Vary S and Time everyday
S <- df$Close
Time <- (end_date - df$Date) / 365

#GBSOption(TypeFlag, S, X, Time, r, b, sigma)@price

df_opt <- rowwise(df) %>%
#this is the premium for one unit of call option  
mutate(premium = GBSOption(TypeFlag = "c",
S = Close,
X = X,
Time = as.numeric((end_date - Date) / 365),
r = r, # interest rate
b = 1.85/100, # dividend yield obtained from https://www.dividend.com/dividend-stocks/financial/investment-brokerage-national/gs-goldman-sachs/
sigma = as.numeric(start_volatility/100))@price,

#this is the delta of a call option (before negation)
delta_hedge = GBSGreeks("delta", TypeFlag = "c", 
                        S = Close, 
                        X = X, 
                        Time =as.numeric((end_date - Date) / 365), 
                        r = r, 
                        b = 1.85/100, 
                        sigma = as.numeric(start_volatility/100))) %>%
ungroup() %>%
  
#delta hedging strategy selected: SHORT CALL LONG STOCK (from BlackS Scholes formula, such strategy should approximate a long position in risk free)
mutate(Option_DoD_PnL = ifelse(Date == start_date, # On the 1st date, we count the cost of buying the option
0, #quantity*premium, #on the first day, receive the call option premium and short the option
-quantity*(premium - Lag(premium))), #if subsequently call option price rises, there is a loss

Hedging_DoD_Pnl = ifelse(Date == start_date, 0, 
                         ifelse(Date == end_date, yes = quantity * Lag(delta_hedge) * (Close - Lag(Close)), 
                                #at the last day, there is no rebalancing of number of shares.
                                no = quantity * delta_hedge * (Close - Lag(Close)))), #long stock - if stock price increase, there is a profit

DoD_PnL = Option_DoD_PnL + Hedging_DoD_Pnl) %>%
mutate(PnL_to_date = cumsum(DoD_PnL),
       HPnL_to_date = cumsum(Hedging_DoD_Pnl), 
       OPnL_to_date = cumsum(Option_DoD_PnL))

maxDrawDown <- {
xs <- df_opt$PnL_to_date
max(cummax(xs) - cummin(xs))
}

# The initial outflow of funds is the cost to buy stocks minus option premium received 
InitialInvt = (df_opt[[1,"delta_hedge"]]*df_opt[[1,"Close"]] - df_opt[[1,"premium"]])*quantity #OUTFLOW of funds
profitability = df_opt[df_opt$Date==end_date,"PnL_to_date"]/InitialInvt
df_opt<-mutate(df_opt, PortValue = InitialInvt + PnL_to_date, PortReturn = DoD_PnL/Lag(PortValue), TTM = end_date- Date)

ggplotly(p=ggplot(df_opt) + geom_line(aes(TTM,Option_DoD_PnL),color = "blue") + ggtitle("option profit - TTM"))
## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.
ggplotly(p=ggplot(df_opt) + geom_line(aes(TTM,Hedging_DoD_Pnl))+ggtitle("stock profit - TTM"))
## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.
kable(tail(df_opt%>%select(PnL_to_date,HPnL_to_date, OPnL_to_date)%>%rename( `final pnl`=PnL_to_date,`hedging pnl`=HPnL_to_date  , `option pnl` = OPnL_to_date ),1))
final pnl hedging pnl option pnl
348.8032 128.7305 220.0727
SR <- as.numeric(((profitability - r/12)/30)/stdev(df_opt$PortReturn, na.rm = TRUE))
cat(paste0("the Sharpe Ratio is ",round(SR,4)))
## the Sharpe Ratio is 0.7669
cat(paste0("The maximum drawdown is ", round(maxDrawDown,4)))
## The maximum drawdown is 348.8032